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We study the dynamic response of a granular chain of particles with a resonant inclusion (i.e., a 
particle attached to a harmonic oscillator, or a mass-with-mass defect). We focus on the response 
of granular chains excited by an impulse, with no static precompression. We find that the presence 
of the harmonic oscillator can be used to tune the transmitted and reflected energy of a mechanical 
, pulse by adjusting the ratio between the harmonic resonator mass and the bead mass. Furthermore, 

we find that this system has the capability of asymptotically trapping energy, a feature that is not 
■ present in granular chains containing other types of defects. Finally, we study the limits of low and 

, high resonator mass, and the structure of the refiected and transmitted pulses. 
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O ■ I- INTRODUCTION 

Granular crystals consisting of tightly packed arrays of solid particles, elastically deforming upon contact, present 
. interesting dynamical features that enabled fundamental physical discoveries and suggested new engineering applica- 

■ tions. Most studies of granular crystals focused on the highly nonlinear dynamic response of these systems |ll-[23|. 
, One-dimensional (ID) sranular crystals have been studied in a number of analytical, numerical and experimental 

• 1-H ■ investigations; see in [2| for a recent review of this topic. The ability to use a wide variety of materials and bead 
sizes, as well as the tunability of the response within the weakly or strongly nonlinear regime makes granular crystals 
a natural paradigm for physical explorations of wave phenomena and the effect of nonlinearity on them [sl, 0] . On 
the engineering side, this tunability makes such crystals promising candidates for numerous applications, including 
shock and energy absorbing materials 1 actuating and focusing devices [1, [ll| , and sound scramblers or filters 

, To model the dynamics of granular chains, Hertzian interparticle interactions (proportional to the relative displace- 

■ ment of adjacent bead centers raised to the 3/2 power in the case of spherical beads) have been established as the 
canonical approach These chains have been shown to support the emergence of nonlinear traveling waves, which 

psj ■ have been described through diff^erent types of partial differential equation models (see, e.g., [ll]) or even by binary 
collision particle models (see, e.g., [l3l)- Although these waves are treated as compactly supported in the continuum 
. approximations, they decay with a doubly exponential power law [23l . [23 | (i.e. extremely fast, but their support is 
' not genuinely compact). 

^ . Once a chain of particles is excited by an impulse, more than 99% of its energy is rapidly and spontaneously 
rearranged into one or more of such nonlinear traveling waves (TW) [l^ . Examining the interaction of the resulting 
TWs with a defect has been of particular interest from the point of view of applications (e.g. for detecting cracks 
or for detecting other irregularities in a medium (see e.g., [26j and references therein). A pioneering study in 
this regard was the computational work of [l^, which examined both the case of a light defect and the resulting 
symmetric emission of traveling waves in both directions by the defect, and that of a heavy defect, which produced a 
train of solitary waves asymmetrically to the right of the defect. The theme of impact upon a light defect has been 
revisited in the work of [l6| where the synergy of experiments, numerical computations and analytical approximations 
demonstrated the possibility of transient breather formation in the system. Very recently, the problem was also 
revisited from a chiefly analytical point of view in [l7| , where a reduction of the problem to a chain of three beads 
solved using a multi-scale expansion was used. This enabled an accurate capturing of the slow dynamics of the defect 
bead and its neighbors, along with the fast transient breathing dynamics of the light bead. It should be noted that 
defects have also been studied in such granular chains in the presence of a precompression force (creating an underlying 
linear limit and hence the potential for a weakly nonlinear regime) in the bulk, where the formation of defect breathing 
modes has been elaborated both analytically/numerically and numerically/experimentally [l9j . The usefulness of 
defect breathing modes (at one edge of the chain) as generators of acoustic diode type effects has also been explored 
recently in [20j . We also mention in passing the consideration of interactions of traveling breathers with defects in a 
"Newton's cradle" model where the Hertzian chain has a local oscillator associated with each bead [27| . 
In this paper, we will consider the interaction of a solitary wave with a new kind of defect, consisting of a secondary 
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mass attached to one of the beads of the chain through a hnear spring, as shown in Fig. [TJ We will refer to this defect 
as a mass with mass (MwM) defect. A MwM defect can be implemented experimentally by attaching a resonant solid 
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FIG. 1: Diagram of a granular chain containing a mass-with-mass defect. 

structure to one of the beads, as shown in Fig. [2l Such a structure should be designed carefully to avoid higher-order 
normal modes of the resonant structure to participate in the dynamics of the system. A ring resonator vibrating in its 
piston normal mode would satisfy this requirement and is able to provide the values of spring stiffness and secondary 
mass needed to reproduce experimentally the effects that we have predicted theoretically. 




FIG. 2: Suggested experimental implementation of the mass-with-mass defect, (a) Defect bead with the ring resonator that 
provides the mass-with-mass. (b) Part of the suggested experimental setup containing the MwM defect. 



Meanwhile, the interaction of the defect bead with the neighboring particles is preserved as Hertzian. A feature of 
a granular chain with a MwM defect is the existence of an underlying linear oscillator at the defect. This oscillator 
presents the potential for long-term energy trapping, a feature that cannot be present in the mass defects considered 
in earlier studies. We examine various limits of the system, especially focusing on the limit of small secondary mass, 
which is analytically tractable and able to reproduce much of the physics of the interaction between a solitary wave 
and a MwM defect. We also separately consider the limit of very large secondary mass, where the reflected pulse has 
much larger energy, as well as the intermediate regime, where an emission of a train of transmitted solitary waves 
that have progressively decreasing amplitudes and speeds is observed. 

Our presentation is organized as follows. In Section |TT1 we introduce the model and the associated quantities that 
will be monitored. In Section IIIIl we discuss the numerical observations of a traveling wave interacting with the 
MwM defect in the regimes of a small, intermediate and large secondary mass. We also briefly comment on the 
case of several defects. In Section llVI we analyze the limit of small secondary mass using two approaches: a direct 
perturbation method that approximates the system in this limit as a local oscillator driven by a solitary wave and a 
multiscale analysis of a reduced model that captures some of the prototypical features of the full system. Section |V] 
contains the summary of our findings and the discussion of possible future directions. 

II. THE MODEL 

We will consider a granular chain of identical spherical elastic beads of mass mi each. Let w„(i) denote the 
displacement of the nth. bead from its equilibrium position and denote Un{t) = u'^^it), Unit) = u'^{t)- The interaction 
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between nth and (n + l)th beads is governed by Hertz interaction potential 



TJ 



x5/2 
Un+l)+ , 



where {x)+ = x when a; > and equals zero otherwise, and a > is a material constant. The beads thus interact only 
when they overlap and are not subject to a force when the overlap is absent. The defect bead, at n = 0, is attached 
to another mass, m2, via a linear spring of stiffness isT > 0, as shown in Fig. [TJ This mass is constrained to move in 
the horizontal direction, with displacement VQ{t). The equations of motion are: 



m2i)o = K{uo - vo), 



3/2 



K{uo - Vo)5nO 



(1) 



where we used the Kronecker delta, (5„o = 1 when n = and zero otherwise. We assume that all masses, except the 
jth bead, for some j < 0, are initially at rest, and the beads in the chain just touch their neighbors at t = 0. The jth 
bead is excited by setting its initial velocity to y > 0. The initial conditions arc thus: 



u„(0) - i;o(0) = 0, 1)0(0) = = u„(0), n^j, uM^V, j < 0. 



(2) 



It is convenient to rescalc ([T|), ([2]) by introducing dimensionlcss displacements w„, vq and time t related to the 
original variables via [l^ 



m2 \ 2/5 



Wo 



T^2\ 2/5 



Wo, 
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The two dimensionless parameters are 



the ratio of the two masses, and 



m2 
mi ' 



K 



^V5^4/5-(/2/5 ' 



(3) 



(4) 



(5) 



which measures the strength of the linear elastic spring in the mass-with-mass defect relative to the Hertzian potential 
stiffness at the particular V . In what follows, we set k = 1 for simplicity (for any we can suitably select V to 
achieve this in Eq. ([5])), while focusing on the effect of e. Substituting ([3]), @ and k = 1 into ([H and © and dropping 
the bars on the new variables, we obtain 



euo = Mo - Wo 



\3/2 



(6) 



and 



M„(0) = wo(0) = 0, wo(0) = = u„(0), n ^ j, 



,(0) = 1, j<0. 



(7) 



We conduct a series of numerical experiments in the spirit of |15l |. in order to understand the dynamics of the 
granular chain in the presence of a mass-with-mass defect. Notice that in addition to the displacement fields u„ and 
vq, and the corresponding velocity fields m„ and vq, another characteristic quantity the system is the total energy 
E = e„, where 



6" = 7;ul + T^'i'l^nO + ^(mo - Vo)'^SnO + ^ 

111 5 



I \5/2 , / ^5/2' 

(Mn-1 - Mn)+ + ("n - ) + 



(8) 



is the energy density (energy stored in each bead). The total energy i? is a conserved quantity of the system. We 
have confirmed that in our dynamical simulations (explicit fourth order Runge-Kutta with fixed time step of 10^'^) 
energy is conserved up to the order of 10^^^. Once the traveling wave interacts with the defect, part of the energy 
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will be reflected, part of the energy will be transmitted, and part of the energy will be trapped. We define the fraction 
of energy that is reflected as 

^ = ;^ E (9) 

n<-l 

the fraction which is transmitted as 

n>l 

so that the trapped fraction of the energy is given by 

l-T-i?=-^e„. (11) 

n=-l 

We now present a detailed discussion of our numerical results and some analytic approximations of the system. 



III. NUMERICAL RESULTS 



We numerically integrate (|6]) for —100 < n < 100 subject to (O, with the defect located at n = 0, in the middle 
of the chain. The initial excitation is at n = j = — 50, which enables the robust formation of the traveling wave well 
before its impact with the mass-with-mass defect site and precludes any backscattering or rebounding waves from the 
excited site to affect the defect location within the duration of our numerical simulations. 

Running the simulations for a range of the mass ratio e = m2/mi values, we compute the fractions of the energy, 
defined in pU|) and (fTT|) . respectively, as functions of e. The results are shown in Fig. [31 A well- understood 




m2/mi m2/mi 

FIG. 3: Left panel: the different fractions of the energy are plotted as a function of the mass ratio e = m2/m\. Right panel: 
trapped portion of the energy. The transmitted energy fraction T of the right part of the chain is presented by the red dash- 
dotted line, the reflected part of the energy R is given by the green dashed line, while the non-vanishing (for e > 0) trapped 
fraction of the energy is shown by the blue solid line. 



limiting case is that of e = 0, when uq = vq, so that the mass-with-mass defect is effectively absent and the system 
features near perfect transmission: T w 1 and i? ~ 0. This observation takes into account the well-known feature [25| 
that over 99% of the impact kinetic energy is stored within a highly nonlinear traveling wave (originally described by 
Nesterenko within his quasicontinuum theory |l|). Another physical limit is that of large values of e, m2 3> mi, when 
the inertia prevents the secondary mass from moving, so ~ 0. In this case the transmission is far from perfect (e.g. 
T ss 0.26 and R ~ 0.71 at e = 10000), and the trapped fraction of the energy approaches the value of about 0.03 as e 
becomes large. Observe that in the regime of comparable bead masses, between these two limits, the trapped energy 
fraction exhibits oscillations. 

We begin by considering the case of small e. As e departs from zero, we still have a single propagating traveling 
wave within the chain, but now it experiences a weak reflection from the defect. We also find that some of the energy 
is trapped, which is a fundamental difference from the observations in the context of light or heavy mass defects in a 
Hertzian chain 
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This picture is corroborated by the detailed dynamics of this case as presented in Fig. |4]for the value of e = 0.1. 
The contour plot in the top left panel presents the space-time evolution of the velocity field which characteristically 
represents the traveling wave. The bottom left panel below the contour plot and on the same time scale shows the 
evolution of the quantities — and (u„ — Un+i) + , which, as we recall, determine the forces exerted on the 

corresponding beads. When the force vanishes, there is what we call a "gap opening" [2l|, i.e., the beads no longer 
interact. It is, thus, clear by also consulting the second and third subplots of Fig. |4] that upon the impact of the 
wave, there emerges a reflection and subsequently a permanent gap opening in the interaction of beads at n = 
and n — — 1. This produces the single reflected pulse present in the process. On the other hand, there is a more 
prolonged interaction between beads at n = and n ~ 1, which, in turn, leads at the early stage of the dynamics to 
the emergence of a transmitted pulse. However, a key feature also arising in the process is the existence of a trapped 
part of the energy. For the small e considered here, this part is weak but it is definitely present and manifests itself as 
oscillations at a characteristic frequency. Interestingly, this frequency, as we will justify in Sec. II VI where we analyze 
this asymptotic case, is precisely the linear frequency of the mass-with-mass oscillator uj = vT+e^. 




FIG. 4: The results of simulations at e = 0.1. The top left panel presents the space-time (n-t) evolution of the contour plot of 
the field u„{i). The thick blue and green lines in the bottom left panel show the quantities (w„_i — u„)+ and {u„ — Un+i) + , 
respectively, to give a sense of when (after the impact) gap openings arise. The top right panel presents the displacements of 
the central sites; cyan dash-dotted line corresponds to u-i{t), thick blue solid to uo{t), green solid to ui{t) and red dashed to 
vo{t). The same pattern is followed for the corresponding velocities at the bottom right plot. 



We now turn to the opposite limit, namely that of a much larger mass TO2 in comparison to the mass mi of the 
granular chain beads. This case is represented by Fig. [Sj which presents the computations for e = 10000. In this case 
the large secondary mass barely moves (dq ~ 0). However, uq oscillates around the zero value after the initial solitary 
wave reaches the defect site, triggering oscillations of ui, whose amplitude descreases with time. As before, the gap 
between n = and n = —I forms, sending a reflected pulse (note, however, the gap is no longer permanent, as these 
two sites briefly interact several times later). As in the small e case, there is also a transmitted pulse in this limit. 
However, the transmitted pulse is much weaker in this case, while the reflected pulse is much stronger: indeed, at 
e = 10000, we have R ~ 0.71 and T « 0.26, far from the near-perfect transmission we observed at small e. It is also 
important to note that there is an immediate gap opening, upon the first passage of the wave, between n — and 
n = 1 in this case, as can be observed in Fig. [5] In this case, the asymptotic scenario (of e — oo) leads the central site 
to acquire a residual oscillation of uq with approximately unit frequency, while vq — >■ 0. In that light, the system to 
consider analytically must include the beads ?i = —1, n — and n = 1 for the u-field only. However, we were unable 
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FIG. 5: The same properties as for the Fig. 3] are illustrated but now for the case of e = 10000. 



to identify an analytical solution of this six-dimensional nonlinear system (note also that unlike the reduced model 
considered in Sec. IIVI for small e, in this case there is no small parameter). Nevertheless, numerical computations 
with the single component system (where there exists a local oscillator at n = 0) clearly confirm the validity of the 
above picture. 

We now examine the interaction of the traveling wave with the defect for intermediate values of e. We start with 
the case e = 10 shown in Fig. |6l when the secondary mass is still considerably larger than the primary one. First, it 
should be noted that as expected here, the fast-scale oscillations are performed by the mass mi (contrary to what was 
the case for small e, where the mass m2 was the one performing the fast-scale vibrations). Second, the phenomena 
observed in this case are significantly different from what we have seen at small e. There still exists a small fraction 
of the energy which remains trapped at the central site. However, the principal phenomenology does not consist of 
a single transmitted and reflected wave, as was the case for small e. In this case, the dominant reflected wave may 
be a single one, as was the case for small e (as shown in Fig. [5]), but there is a cascade of transmitted waves which 
are somewhat reminiscent of the corresponding train of transmitted solitary waves discussed in 15| for large mass 
defects in Hertzian chains; see also the work of [l^] for a similar problem where a solitary wave train also emerges 
as a consequence of a chain being stricken by a heavier bead, a role that in our system is played by the defect bead. 
The top right panel of the figure illustrates the displacements and allows the observation of the gap openings where 
there is no force after the impact of the traveling wave. We see that upon the transmission of the original traveling 
wave, there is a gap opening between the defect site and its left neighbor (n = and n = — 1; see the thick blue 
curve in the bottom left panel in Fig. |6]), which initiates the principal reflected wave traveling to the left. These two 
sites never interact again. However, the interaction of the n = defect site with the site to the right of it (n = 1) is 
different. What is observed is that there is a series of gap openings which are followed (each in turn) by subsequent 
compression intervals (see the thin green curve in the bottom left panel). Each one of these sequences produces a 
new traveling wave, as can be seen in the top left panel. However, naturally, as the energy trapped within the central 
site keeps decreasing upon the release of subsequent traveling waves, the amplitude of each later (emerging) wave in 
the sequence is weaker than that of its predecessors and hence its speed is also smaller. This, in turn, justifies this 
sequence of decreasing amplitude and speed traveling waves emitted by the defect. 

We note that a similar sequence of gap openings between the beads at n = and n = 1 was also seen at e = 0.1 (see 
the thin green curve in the inset in the bottom left panel Fig. |4]). However, due to the much smaller amplitude and 
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higher frequency of oscillations of uo(t) in that case, the nonzero interaction forces during the compression intervals 
and the duration of these intervals were too small to initiate additional traveling waves. 
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FIG. 6: The same properties as for the Fig. |4] are illustrated but now for the case of e = 10. Here, it is the lighter mass mi 
executing the fast scale oscillations, while opening a gap from the site to its left and presenting a cascade of alternating gap 
openings and compression events with its right neighbor which, in turn, lead to the successive emission of progressively smaller 
amplitude (and thus slower) waves within the solitary wave train that emerges in the space-time contour plot. 

The dynamics for the case when the two masses are comparable is represented in Fig. [71 which corresponds to the 
value 6 = 0.5. It should be noted here that this case represents the parameter close to the value at which the trapping 
fraction of the energy is maximal. In particular, as observed in Fig. [3l there is a clearly defined maximum in the 
fraction of energy that can be trapped by the system; this is another unique feature of our system. Not only can 
energy be trapped by the mass-with-mass defect but the trapping has a non-monotonic dependence on the ratio of 
the masses. As illustrated in the figure, here the scenario is different than the ones observed before. Here, the defect 
bead does not lead to a gap opening with respect to the bead to its left, but rather (predominantly) with respect 
to the bead to its right. In this case, we observe a two-peaked structure in (ui — mo)+ (see the bottom left panel of 
Fig. [7]) which appears to nucleate two traveling waves moving in the right half of the chain. While a gap indeed opens 
past this principal interaction between n = and n = 1, the beads do interact anew at a much later time (around 
t = 110). Yet, similarly to the interactions shown in Fig. lU such later (and considerably weaker) interactions do not 
produce appreciable solitary waves traveling to the right of n = 0. 

Finally, we turn to a setting where there exist additional defects within the chain. This is shown in Fig. [51 which 
contains both the case of two adjacent MwM defects (left panel), as well as that with three such adjacent defects (right 
panel). It can be clearly discerned that the addition of further defects significantly enhances the trapped fraction of 
the original energy within the defect sites. This fraction increases by about 10% (for comparable values of 111.2/ mi), 
when the defective region expands from one to two sites, and a comparable increase is observed when a third site is 
added. It is observed that this increase is at the expense of the transmitted fraction, which is well below 20% (« 15% 
and well below the trapped fraction) for the three-MwM defect case. Interestingly, the trapped fraction preserves its 
oscillatory structure (notice that a less pronounced, yet somewhat similar variation is observed in the transmitted 
fraction), but the peak locations vary, as the number of defect sites is increased; e.g. the first, most pronounced peak 
occurs at larger values of 1112/1711. 





FIG. 8: The same diagnostics as those of Fig. |3] are shown here (transmitted T, reflected R and trapped 1 —T — R fractions of 
the energy) as a function of the mass ratio, but now for the case of two adjacent mass defects (left) and that of three adjacent 
mass defects (right). 



IV. ANALYSIS OF THE SMALL MASS RATIO CASE 



We now present a theoretical formulation which enables a semi-quantitative understanding of the situation at hand 
for the case of small e. We will use two different approaches, one of which involves perturbation analysis of the full 
system and the other is based on multiscalc analysis of a reduced system. 
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A. Perturbation analysis 

Assume that there exists a traveling wave satisfying the standard Hertzian chain equation u„ — (w„-i — Un)V^ 



{un — u„+i)ij/^, where the solution can be well approximated by the Nesterenko solitary wave [l[ with 




2 \ n 

-(n-c(t-to)) j, to + ^ <t<to + -, (12) 

otherwise 

where — (4/5)^^/^ and to is the time when the wave leaves the n, = site. Now, focusing on the difference 
Wq = uq — vq, we note that it satisfies the equation 

Wo + uj'^wo = (w_i - Mo)+ ^ - [uo - ^, (13) 

where 

is the natural frequency of the MwM defect. While Eg.f [T^ is exact, observe that to a leading order its right hand 
side can be approximated by uq. This means that the linear oscillator is driven by the propagating traveling wave, 
which we further approximate by the above Nesterenko expression (|12p at 7i = 0, with uq obtained by differentiating 
([T^ once with respect to time. Notice that this drive is only active within the time interval [to — (7r-\/5/2/c), to]: i-G- 
as the wave passes over n = site. 

Within this approximation, Eq. ([T3|) can be solved exactly to give 

/■* sin(i:j(t - r)) .. 

Wo = / -Uo{T)dT. (14) 

Jo ^ 

This simplified model has the advantage that it can be evaluated explicitly: 

'0, i<^o-^ 

(^^/T0(5w2 „ 32^2) sin(20) + -y/KSc^ - 5^^^) ^-^^^^^ _^ sin(t^(^ + t - to})^ , to - ^ < t < to 

_i9gdcisin(a.i^)cos(^(t-to + ^)), t > to, 

(15) 

where </> = \/2/5(c(t — to)) and D ~ 256cf^ — 200c^lo^ + 25uj^. Note that it corroborates the numerical observation 
that the trapped oscillation of the mass-with-mass defect is executed with the natural frequency lo of the relevant 
oscillator. 

The results of this approximation for small values of e are illustrated in the example of Fig. ^ Here the traveling 
wave is initiated by the standard initial condition above and its resulting speed c and amplitude A are computed from 
the simulation prior to its impact with the mass-with-mass defect (we obtain c ~ 0.841 and A « 0.811). Then the 
time derivative of the expression of Eq. (jisp (shown by the green line) is evaluated and compared, after an appropriate 
time shift to, with the direct result of the numerical simulation (black line) at the central site. It can be seen that 
the analytical solution captures the numerical result quite well qualitatively and, as expected, properly captures the 
frequency of the residual vibration trapped at the central site. It does not, however, capture its amplitude in general, 
which can be either larger (see, for example, the discrepancy at e = 0.005) or smaller (e.g. e = 0.025) than the 
numerical value. It also underestimates the amplitude of the large pulse that precedes it. This happens because our 
approximation neglects the effect the defect has on the propagating wave which adjusts its speed and shape prior to 
reaching the n = site. 



B. Small e limit: two-bead problem 



To better understand the small e limit, we now turn to a simplification of the original problem in which we only 
consider the system involving two beads. Neglecting the left part of the chain can be justified by the fact that the 
numerical simulations show that after the defect bead has been kicked by the previous one, there is a gap opening 
between the two and therefore they do not interact again. Hence, the simplest configuration that could capture the 
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transmission of the wave from n = to n = 1, as well as the trapping of the energy in the n = site would be the 
two-site setting examined below. Importantly, these two phenomena exhibit a separation of time scales: the oscillation 
within the defect bead is one of a fast time scale, while the interaction between n = and n = 1 occurs in a slower 
time scale, enabling the multiscale analysis presented below. 
The simplified dynamical equations then read: 



ilo = -{uq - uifj^ - {uq - Vq) 

I \3/2 



(16) 



In considering this reduced model, we follow the approach in [17[ for ^ problem with a light mass defect, with the 
goal to qualitatively capture the dynamics of the system during an initial time period after the propagating wave 
reaches the mass-with-mass defect. For example, consider the initial conditions (emulating the kick provided through 
the gap-opening interaction with rt = — 1) 

uo(0)=i-o(0)-ui(0)-i)o(0)-?ii(0) = 0, uo(0) = l. (17) 

The mmierical solution of the system p6)) at e = 0.01 and e = 0.1 subject to (|17p is shown in Fig. [TOl By numerically 



e=0.01 




6=0.1 




FIG. 10: Numerical solution of (HH), ((TZl) at e = 0.01 and e = 0.1. 



integrating the simplified system and comparing the results with the previously obtained ones for the full system, we 
see that the simplified model qualitatively captures the behavior of the full system at small e in that vq (t) and uq (t) 
both oscillate, vo{t) with a larger amplitude, about some average value that increases with time, initially rapidly and 
then more slowly. 
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To analytically approximate p6|) at small e, we use a two-timing approach [T^]- We introduce the fast time r = i/v^. 
This rescaling of time yields: 



d^UQ 3/2 , , 

-j-^ = -e{uo - ui)^ - e{uo - vq) 



Uq - Vo 



dr2 



(18) 



dr2 

Substituting the expansion 

Mr) = Vo^°^ (r, t) + V~eV^'^ (r, t) + eV,^''> (r, t) + 6^/ Vj'^ (r, + 0(e2) 
ui(r) = u[°\T,t) + V~eU['\T,t)+eU[^\T,t)+€^/^u[^\T,t) + 0{e^) 

in (jl8p . expanding the right hand side in terms of e and keeping the terms up to O(e^), we obtain 



Ot^ ^ \ dr^ drdt I \ dt^ dtdr dr^ I \ dt^ 









drdt J 












drdt J 







r(dH^ ^^v^\ (^v^ ^^v^ ^V^\ yjd^ 

dT^ + 1^ + drdt J^'y dt^ ^ dtdr + 9r2 j + ' dt^ 









dtdr 












dtdr 


9r2 


a2[/f) 

dtdr 


a2[/f) 

^ 9r2 



Solving the above equations and eliminating the secular terms (which feature unbounded growth in r), we obtain 

U^'^^Boit), ui"^=B,{t) (20) 

and a fast oscillation of the defect as 

^(1°) = C{t) COST + D{t) sinr + Bo{t). (21) 

Next we consider the 0{y/e) problem: 

32rr(l) ^,2rr(0) a2rr(l) a2rr{0) ;.2t/(1) n2T/(0) 



^ ^0 12 ° -0 ^ ^1 12 ^ -0 ^ 12 ° - L^<'^ V^^^ (22) 



From it follows that 



92^(°) _ 92^(0) 
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and hence the first two equations in (1221) imply that 

C/W=i?oW, ui'^=E,{t). (23) 
Differentiating (PT|) to find the mixed derivative in the third equation in (P^. 

g2y{0) 

" =- C" (t) sin T + i?'(t) COST, 



drdt 



we find that we must have C'(t) = D'{t) = in order to avoid seeular terms in Vq^\ Thus C and D are constant, 
and we have 

Vo^°^ = C cos T + L» sin t + Bo(t), (24) 
while the last equation in ([^^ and ([231) yield 

^((1) = i^(t)cosT + G(t)sinT + £;o(0- (25) 
In the subsequent order, namely 0(e), we have 





f 2 ° 






f 2 ° - 


9r2 

a2y(^) 


a2C/f) 


dtdr 
f 2 ^ 


a2[/f) 


9*2 




9r2 



+ 

Using ((20|) . (pS]) and recalling the third equation in (fT9)) . we can rewrite the first equation in (j26p as 



,2 



9x2 



+ V^'^) = -B^{t) {Bo{t) B,{t)ff. 



(2) 

To eliminate secular terms from Uq , we must then set the right hand side of this equation to zero: 



B'^{t)^-iBoit)-B,it)ff. (27) 



Elimination of the linear term in r then means is a function of t only, yielding 

r(2) 



U^'=Hoit)-CcosT-DsmT, (28) 



ul"^' , we obtain 



where we used ([24|) . Similarly, using ([201) and psp in the last equation of ([26]) and eliminating the secular terms in 

B';it)^iBoit)-B,it)ff (29) 



and 

C/p^=ffi(t). (30) 
Considering now the second equation in ([26]) and substituting the mixed derivative of (|25p and ([28]) . we have 



+ ^ _(2G'(t) + c) COST + (2F'(t) - D) sinr + i7o(t) - B'^it). 

To avoid secular terms, we must set the coefficients in front of sinr and cost in the right hand side to zero, which 
yields F{t) = {D/2)t + Dq and G{t) = —{C/2)t + Cq, where Co and Z?o are constant. Thus we have 



^o'^ ^ [ + Do] COST + [Ca - ^t] siuT + Eo[t), (31) 
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while 



Vo(2) ^ j(^) COST + J{t) sinr + Ho{t) - B'^(t). 



Finally, wc turn to the O (e'^^'^) problem: 



9t2 


f 2 - 


ar2 




f 2 " - 








(9r2 




f 2— ^ - 


92[/p) 

1- ^ „ 



2 

C/(^)-K<^) (32) 



9t2 dtdr dr^ 2 



(1) 



Proceeding as in the 0{e) case, we use (f^ . ([25]), the third equation in ([^ and ([M)) to rewrite the first equation in 
(1321) as 



5 ,,fn\ ^/,/,^ 3 

9^2 

("3) 

To eliminate secular terms from Uq , we must then set the right hand side of this equation to zero: 



Kit) = (Boit) - B,{t))Y' (Eoit) - E,it)). (33) 
Then, again eliminating the secular terms, we get that Uq^"^ + Vq^"^ is only a function of t, which in light of ([5T|) yields 

U^^'^ = Ko{t) - (jt + Do) COST - (co - ^t) siuT. 



Similarly, ([23]), ([30]) and the last equation in (|32]) result in 



EUt) = I (Boit) - i?i(t))f (i?o(i) - E,{t)) (34) 



and C/f ^ = Ki{t). 

So up to order 0(e) we have 



Co - ) sin —= )■ , 



Mo(t) - Bo{t) + VeEoit), ui{t) = Bi{t) + ^feEl{t), 
vo{t) = Bo{t) + C cos ^ + D sin + y/I I Eo{t) + ( Dq + ^t] cos ^ 

where Bo{t), Bi{t) are found by solving ([27]), ([^ . Ef){t), Ei{t) satisfy ([55]) . ([M]) . and the initial conditions for these 
hmctions, as well as the constants Dq, D, Cq, C, are found from the initial conditions for uo{t), ui{t) and va{t). 

Now we consider the initial conditions (|T7]) . Then uo(0) = vo{0) = ui(0) = 0, which imply that i?o(0) = Bi{Q) = 
Eo{0) = Ei{0) = and C = L»o = 0. The initial condition uo(0) 1 implies B^(0) = 1 and £'^(0) = 0, while 
1*1(0) = implies B[{0) ^ E'i{0) = 0. Differentiating vo{t) and using uo(0) = and the above results, we get 

^o(O) = l + -^ + Co + ^=0, 
ye Z 

which implies D = and Cq — —1. Since EQ{t), Ei{t) satisfy the linear system (|33p . (|34| subject to the zero initial 
conditions, they must vanish. Thus we have 

uo{t)^Bo{t), ui{t)^B,{t), Vn{t)^Bo{t)-^esm^ (35) 

where Bo{t), Bi{t) satisfy ([17]), ([H]) and the initial conditions Bo(0) = Bi(0) = Bj(0) = and 5^,(0) = 1. To solve 
this problem, it is convenient to introduce the new variables yo(0 = {Boif) — Bi{t))/2 and zo(i) = {Bo{t) + Bi[t))/2. 
Then zo(t) satisfies ZQ{t) = 0, zo(0) = 0, Zq{Q) = 1/2, so we have z{t) = t/2. Meanwhile, 2/o(t) solves 

2/^,' + (2yo)+/' = 0, 2/0(0) = 0, 2/^,(0) = 1/2. (36) 
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When yo > 0, it satisfies j/q + (2j/o)'^^^ = 0, which together with the initial conditions imphes that it hes on the 
trajectory 

This can be solved for t as a function of in terms of hypergeometric functions. We obtain 



t = 



Q<t<U 



2y,,F,[lllfV2yl'\ 



(37) 



where 



22/0 2F1 



[2 17 

1-5' 2' 5 



, 1] is such that yo(^*) 



32^/2 



2/5 



?/^(2t^) = -1/2 and j/o(2i*) = 0, this yields 

yo{t)=U-t/2 ioTt>2U. 
Recalling that Bo{t) = yo{t) + zo{t) and Bi{t) = zo{t) — yo{t), we obtain 

Uo{t) 



For t > 2t^ we have j/g = 0, and since 

(38) 



i+yo{t)i ■Ui{t)^^-yo{t), vo{t) = ^ +yo{t) - ^/esin^ 
2 2 2 ^/6 



(39) 



where yo{t) is given by p7p . (j38p . In particular, for t > 2t^ we have uq becomes constant, uo(i) = i*, while ui linearly 
increases, ui{t) = + i/2. 

The comparison of the approximation p9p (black curves) and the numerical solution (colored curves) is shown in 
Fig. [TT]at e = 0.01. One can see that the approximation works very well for t < 5, but starts deviating for larger 

e=0.01 



rt/"' "o 




FIG. 11: Comparison of the numerical solution of (|16|) . (|17p at e = 0.01 (colored solid curves) and the two-timing approximation 
up to 0(e) (black dashed curves). 



t because it does not capture the slight increase in uo{t) for t > 2t,. It also does not capture the small-amplitude 
oscillations of ua{t), which become visible at larger e; see, for example. Fig. [10] at e = 0.1. To capture these effects, 

one needs to include higher-order terms. Note, however, that the function (uq — wi)+ ^ is not smooth at zero and 
cannot be expanded beyond the first derivative term at this point. Hence, while this approach is valuable in analytical 
understanding of the leading order dynamics in the simplified two-site system, it also has its limitations with respect 
to some of the higher order effects therein. 

Nevertheless, the analytical approximations considered in this section offer a detailed quantitative picture of the 
energy trapping and residual oscillatory dynamics within the defect (and capture its frequency), as well as of the 
detailed exchange dynamics between sites n = and n = 1. 



V. CONCLUSIONS AND FUTURE CHALLENGES 



In the present work, we have explored the propagation of waves in a granular chain in which a mass-with-mass 
defect is present. In this setting, a local oscillator with additional parameters arises, the most significant of which is 
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the ratio between the mass-with-mass defect and that of the rest of the masses within the chain. We have considered 
the problem as a function of this parameter numerically and wherever possible also analytically. 

We found that in the case of the small defect to bead mass ratio, the traveling wave remains essentially unaltered, 
except for the fact that a part of its energy is reflected and a part of its energy is trapped in the form of localized 
oscillation. The trapping of energy is an exclusive feature of this system that was not observed in chains contain- 
ing mass defects. This phenomenon was studied analytically in two distinct ways. One of them involved a direct 
perturbative approach based on the fact that the local oscillator is principally driven by the weakly affected (in this 
case) solitary wave. The second one was a multiscale technique applied on a reduced, two-beads system based on 
two-timing which revealed the key role of the dynamics/interaction of the defect site n = with the following one 
(n = 1). We also studied the effect of the defect for all the representative values of the mass ratio using numerical 
methods. In the limit of very large mass ratio, we found that the reflection is more significant than the transmission 
and a considerable amount of trapping still occurs. In the intermediate mass ratio cases we found the potential of 
exciting multiple waves either one directly after the other, or through a sequence of gap openings leading to a train of 
solitary waves emitted towards the right of the defect. In each case, we studied the dynamics of the interaction with 
the defect by computing the trajectories of the beads in the vicinity of the defect (n = — 1, n = and n = 1). 

The present study suggests many possible future investigations on systems containing MwM defects. On one hand, 
the large mass ratio limit would certainly be worthwhile to consider analytically. We also found that the energy 
trapped in the MwM defect shows a non-monotonic dependence on the mass ratio which has not yet been explained. 
And finally, this study could be extended to cover the propagation of traveling waves in ordered systems where all the 
beads have a MwM defect. These topics are currently under consideration and will be reported in future publications. 
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